{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Integrating neural models using exact integration \n",
    "\n",
    "\n",
    "## The simple integrate-and fire model\n",
    "\n",
    "For the simple integrate-and-fire model the voltage $V$ is given as a solution of the equation:\n",
    "\n",
    "$C\\frac{dV}{dt}=I$.\n",
    "\n",
    "This is just the derivate of the law of capacitance $Q=CV$. When an input current is applied, the membrane voltage increases with time until it reaches a constant threshold $V_{\\text{th}}$, at which point a delta function spike occurs.\n",
    "\n",
    "A shortcoming of the simple integrate-and-fire model is that it implements no time-dependent memory. If the model receives a below-threshold signal at some time, it will retain that voltage boost until it fires again. This characteristic is not in line with observed neuronal behavior.\n",
    "\n",
    "## The leaky integrate-and fire model\n",
    "\n",
    "In the leaky integrate-and-fire model, the memory problem is solved by adding a \"leak\" term $\\frac{-1}{R}V$ ($R$ is the resistance and $\\tau=RC$) to the membrane potential:\n",
    "\n",
    "(1)   $\\frac{dV}{dt}=\\frac{-1}{\\tau}V+\\frac{1}{C}I$.\n",
    "\n",
    "This reflects the diffusion of ions that occurs through the membrane when some equilibrium is not reached in the cell.\n",
    "\n",
    "\n",
    "## Solving a  homogeneous linear differential equation\n",
    "\n",
    "To solve (1) we start by looking at a simpler differential equation:\n",
    "\n",
    "$\\frac{df}{dt}=af$, where $f:\\mathbb{R}\\to\\mathbb{R}$ and $a\\in\\mathbb{R}$\n",
    "\n",
    "Here the solution is given by $f(t)=e^{at}$.\n",
    "\n",
    "## Solving a non-homogeneous linear differential equation\n",
    "When you add another function $g$ to the right hand side of our linear differential equation:\n",
    "\n",
    "$\\frac{df}{dt}=af+g$ (this is now a non-homogeneous differential equation)\n",
    "\n",
    "things (can) become more complicated.\n",
    "\n",
    "### Solving it with variation of contants\n",
    "\n",
    "This kind of differential equation is usually solved with \"variation of constants\" which gives us the following solution:\n",
    "\n",
    "$f(t)=e^{ct}\\int_{0}^t g(s)e^{-cs}ds$.\n",
    "\n",
    "This is obviously not a particularly handy solution since calculating the integral in every step is very costly.\n",
    "\n",
    "### Solving it with exact integration\n",
    "\n",
    "With exact integration these costly computations can be avoided. \n",
    "\n",
    "#### Restrictions to $g$\n",
    "But only for certain functions $g$! I.e. if $g$ satisfies (is a solution of):\n",
    "\n",
    "$\\left(\\frac{d}{dt}\\right)^n g= \\sum_{i=1}^{n}a_i\\left(\\frac{d}{dt}\\right)^{i-1} g$ \n",
    "\n",
    "for some $n\\in \\mathbb{N}$ and a sequence $(a_i)_{i\\in\\mathbb{N}}\\subset \\mathbb{R}$.\n",
    "\n",
    "For example this would be the case for $g=\\frac{e}{\\tau_{syn}}t e^{-t/\\tau_{\\text{syn}}}$ (an alpha funciton), where $\\tau_{\\text{syn}}$ is the rise time.\n",
    "\n",
    "### Reformulating the problem\n",
    "\n",
    "The non-homogeneous differential equation is reformulated as a multidimensional homogeneous linear differential equation:\n",
    "\n",
    "$\\frac{d}{dt}y=Ay$ where \n",
    "\n",
    "$A=\\begin{pmatrix}\n",
    "a_{n}&a_{n-1}&\\cdots&\\cdots&a_1&0\\\\\n",
    "1&0&\\cdots&0&0&0\\\\\n",
    "0&\\ddots&\\ddots&\\vdots&\\vdots&\\vdots\\\\\n",
    "\\vdots&\\ddots&\\ddots&0&0&0\\\\\n",
    "0&0&\\ddots&1&0&0\\\\\n",
    "0&0&\\cdots&0&\\frac{1}{C}&-\\frac{1}{\\tau}\\\\\n",
    "\\end{pmatrix}$\n",
    "\n",
    "by choosing $y_1,...,y_n$ canonically as:\n",
    "\n",
    "$\\begin{align*}\n",
    "y_1&=\\left(\\frac{d}{dt}\\right)^{n-1}g\\\\\n",
    "\\vdots&=\\vdots\\\\\n",
    "y_{n-1}&=\\frac{d}{dt}g\\\\\n",
    "y_{n}&=g\\\\\n",
    "y_{n+1}&=f.\n",
    "\\end{align*}$\n",
    "\n",
    "This makes ist very easy to determine the solution as\n",
    "\n",
    "$y(t)= e^{At}y_0$ and \n",
    "\n",
    "$y_{t+h}=y(t+h)=e^{A(t+h)}\\cdot y_0=e^{Ah}\\cdot e^{At}\\cdot y_0=e^{Ah}\\cdot y_t$.\n",
    "\n",
    "This means that once we have calculated $A$, propagation consists of multiplications only.\n",
    "\n",
    "### Example: The leaky integrate and fire model with alpha-funcition shaped inputs (iaf_psc_alpha)\n",
    "\n",
    "The dynamics of the membrane potential $V$ is given by:\n",
    "\n",
    "$\\frac{dV}{dt}=\\frac{-1}{\\tau}V+\\frac{1}{C}I$\n",
    "\n",
    "where $\\tau$ is the membrane time constant and $C$ is the capacitance. $I$ is the sum of the synaptic currents and any external input:\n",
    "\n",
    "Postsynaptic currents are alpha-shaped, i.e. the time course of the synaptic current $\\iota$ due to one incoming spike is\n",
    "\n",
    "$\\iota (t)= \\frac{e}{\\tau_{syn}}t e^{-t/\\tau_{\\text{syn}}}$.\n",
    "\n",
    "The total input $I$ to the neuron at a certain time $t$ is the sum of all incoming spikes at all grid points in time $t_i\\le t$ plus an additional piecewise constant external input $I_{\\text{ext}}$:\n",
    "\n",
    "$I(t)=\\sum_{i\\in\\mathbb{N}, t_i\\le t }\\sum_{k\\in S_{t_i}}\\hat{\\iota}_k \\frac{e}{\\tau_{\\text{syn}}}(t-t_i) e^{-(t-t_i)/\\tau_{\\text{syn}}}+I_{\\text{ext}}$\n",
    "\n",
    "$S_t$ is the set of indices that deliver a spike to the neuron at time $t$, $\\tau_{\\text{syn}}$ is the rise time and $\\iota_k$ represents the \"weight\" of synapse $k$.\n",
    "\n",
    "#### Exact integration for the iaf_psc_alpha model\n",
    "\n",
    "First we make the substitutions:\n",
    "\n",
    "$y_1=\\frac{d}{dt}\\iota+\\frac{1}{\\tau_{syn}}\\iota$\n",
    "\n",
    "$y_2=\\iota$\n",
    "\n",
    "$y_3=V$\n",
    "\n",
    "for the equation \n",
    "\n",
    "$\\frac{dV}{dt}=\\frac{-1}{Tau}V+\\frac{1}{C}\\iota$\n",
    "\n",
    "we get the homogeneous differential equation (for $y=(y_1,y_2,y_3)^t$)\n",
    "\n",
    "$\\frac{d}{dt}y= Ay=\n",
    "\\begin{pmatrix}\n",
    "\\frac{1}{\\tau_{syn}}& 0 & 0\\\\ \n",
    "1 & \\frac{1}{\\tau_{syn}} & 0\\\\ \n",
    "0 & \\frac{1}{C} & -\\frac {1}{\\tau}\n",
    "\\end{pmatrix}\n",
    "y$.\n",
    "\n",
    "The solution of this differential equation is given by $y(t)=e^{At}y(0)$ and can be solved stepwise for a fixed time step $h$:\n",
    "\n",
    "$y_{t+h}=y(t+h)=e^{A(t+h)}y(0)=e^{Ah}e^{At}y(0)=e^{Ah}y(t)=e^{Ah}y_t$.\n",
    "\n",
    "The complete update for the neuron can be written as\n",
    "\n",
    "$y_{t+h}=e^{Ah}y_t + x_{t+h}$\n",
    "\n",
    "where \n",
    "\n",
    "$x_{t+h}+\\begin{pmatrix}\\frac{e}{\\tau_{\\text{syn}}}\\\\0\\\\0\\end{pmatrix}\\sum_{k\\in S_{t+h}}\\hat{\\iota}_k$\n",
    "\n",
    "as the linearity of the system permits the initial conditions for all spikes arriving at a given grid point to be lumped together in the term $x_{t+h}$. $S_{t+h}$ is the set of indices $k\\in 1,....,K$ of synapses that deliver a spike to the neuron at time $t+h$.\n",
    "\n",
    "The matrix $e^{Ah}$ in the C++ implementation of the model in NEST is constructed [here](https://github.com/nest/nest-simulator/blob/b3fc263e073f46f0732c10efb34fcc90f3b6771c/models/iaf_psc_alpha.cpp#L243).\n",
    "\n",
    "Every matrix entry is calculated twice. For inhibitory post synaptic inputs (with a time constant $\\tau_{syn_{in}}$) and excitatory post synaptic inputs (with a time constant $\\tau_{syn_{ex}}$).\n",
    "\n",
    "And the update is performed [here](https://github.com/nest/nest-simulator/blob/b3fc263e073f46f0732c10efb34fcc90f3b6771c/models/iaf_psc_alpha.cpp#L305). The first multiplication evolves the external input. The others are the multiplication of the matrix $e^{Ah}$ with $y$. (For inhibitory and excitatory inputs)\n",
    "\n",
    "\n",
    "\n",
    "## References:\n",
    "\n",
    "[1] RotterV S & Diesmann M (1999) Exact simulation of time-invariant linear\n",
    "    systems with applications to neuronal modeling. Biologial Cybernetics\n",
    "    81:381-402; which you will find [here](http://dx.doi.org/10.1007/s004220050570)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
